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Abstract 


This  report  provides  a  summary  of  the  work  carried  out  in  the  second  year  of  the  SOARD  project, 
Grant  No.  FA9550-15-1-0069,  devoted  to  the  investigation  and  improvement  of  the  detection  and 
tracking  methods  of  inactive  Resident  Space  Objects  (RSOs). 

In  the  second  year,  a  Random  Finite  Set  (RFS)  based  Joint  Target  Detection  and  Tracking  filter 
was  evaluated  for  the  space  object  tracking  scenarios  and  two  extensions  were  developed  in  order  to 
increase  their  robustness  to  unknown  detection  statistics.  The  performance  of  the  extensions  was 
evaluated  using  both  simulated  and  real  data  from  the  Chilbolton  Advanced  Meteorological  Radar 
(CAMRa).  Both  pre-processed  and  raw  data  sets  from  the  CAMRa  were  obtained.  In  the  pre- 
processed  data  set,  a  maximum  of  one  detection  was  reported  per  bearing  angle  based  on  detection 
methods  used  at  the  CAMRa  site.  The  raw  data  was  also  processed  by  using  a  Constant  False 
Alarm  Rate  algorithm  that  was  able  to  report  multiple  detections  per  bearing  angle,  increasing  the 
probability  of  detecting  the  target  of  interest.  Tracking  algorithms,  based  on  the  RFS-based  Joint 
Target  Detection  and  Tracking  (JoTT)  filter  were  investigated,  which  also  include  the  estimation  of 
target  detection  statistics,  parameters  which  are  often  difficult  to  determine,  and  which  can  be  time 
varying.  In  all  cases,  qualitative  and  quantitative  improvements  in  the  robustness  of  the  proposed 
tracking  methods  were  observed  when  comparing  with  the  standard  JoTT  filter. 

In  addition,  image  sequences  from  the  Georgia  Tech  observatory,  known  as  the  Omnidirectional 
Space  Situational  Awareness  (OmniSSA)  data  set,  were  processed  to  determine  the  feasibility  of 
applying  the  RFS  filtering  concepts  to  image  data. 
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1  Introduction 


The  growing  number  of  inactive  RSOs  poses  a  real  threat  to  future  space  missions  and  has  provoked  an 
increased  need  for  developing  robust  methods  of  orbital  object  tracking.  As  an  example,  the  Chinese 
anti-satellite  weapon  test  in  2007  and  the  unintentional  collision  of  the  Cosmos-2251  and  Iridium-33 
satellites  in  2009  are  estimated  to  have  produced  almost  36%  of  the  total  number  of  currently  tracked 
objects  in  the  Low-Earth  Orbit  (LEO)  [18].  Currently,  the  number  of  objects  being  tracked  by  the  U.S. 
Space  Surveillance  Network  (SSN)  exceeds  21,000  targets  and  mainly  consists  of  objects  which  have  an 
effective  diameter  larger  than  10cm  [18]. 

Detecting  and  tracking  these  objects  is  of  utmost  importance  to  avoid  potential  collisions  for  future 
space  missions.  The  following  challenges  arise  in  the  orbital  object  tracking  problem: 

•  The  highly  non-linear  nature  of  the  motion  models  and  the  presence  of  a  number  of  external 
forces,  other  than  gravitational  interaction,  such  as  solar  drag  and  third-body  perturbations,  adds 
computational  complexity  to  the  state  propagation  models. 

•  Problem-specific  Initial  Orbit  Determination  (IOD)  is  necessary  when  the  measurements  of  the  full 
state  (i.e.  both  position  and  velocity)  of  the  object  required  for  propagation  are  not  available. 

•  Due  to  the  relatively  small  observation  times  and  the  limited  field  of  view  of  a  single  sensor,  the 
object  remains  unobserved  most  of  the  time. 

•  The  problem  of  noise  and  clutter  is  common  to  all  tracking  problems,  however,  combined  with 
the  previously  mentioned  issues,  tracking  RSOs  poses  more  challenges  than  traditional  tracking 
scenarios. 

•  The  relatively  low  illumination  of  debris  compared  to  other  objects  such  as  stars  or  the  moon  further 
complicates  their  detection  and  code. 

In  this  report,  the  work  carried  out  during  this  final  period  is  shown.  First,  the  robust  versions  of 
the  filter  implemented  in  year  1  for  radar  data  are  shown.  These  extensions  do  not  rely  on  particular 
values  for  the  parameters  of  the  detection  statistics.  Instead,  these  parameters  are  incorporated  into  the 
filtering  process.  Second,  initial  work  towards  the  use  of  visual  data  is  presented.  To  this  end,  a  data  set 
obtained  at  Georgia  Institute  of  Technology,  the  OmniSSA  data  set,  is  used.  Initial  work  using  this  data 
set  corresponds  to  the  required  pre-processing  of  images,  detection  of  astronomical  objects  such  as  stars 
in  order  to  obtaining  the  corresponding  mapping  between  image  and  celestial  coordinates. 

This  report  is  organized  as  follows.  In  Section  3  a  description  of  the  radar  data  set  used  to  evaluate 
the  proposed  robust  filtering  algorithms  and  the  visual  data  set  used  to  evaluate  detection  methods  are 
shown.  Section  4  presents  the  robust  extensions  to  the  standard  JoTT  filter  in  order  to  be  less  dependent 
on  particular  values  of  the  detection  statistics’  parameters.  Then,  in  Section  5  the  initial  procedures  are 
carried  out  using  the  visual  data  set  are  described.  Section  6  shows  the  results  of  the  newly  proposed 
filtering  algorithms  and  of  the  initial  procedures  performed  in  the  visual  data,  and  finally  in  Section  7 
possible  future  work  is  discussed. 

2  Project  Outreach 

This  project  has  resulted  in  two  published  articles,  one  conditionally  accepted  publication,  and  one 
workshop  presentation  as  follows: 

•  ISI  Journal  Article:  v Robust  Joint  Target  Detection  and  Tracking  for  Space  Situational  Awareness ”, 
Andrey  Pak,  Javier  Correa  and  Martin  Adams,  AIAA  Journal  of  Guidance,  Control,  and  Dynamics, 
Conditionally  Accepted  for  publication. 

•  ISI  Journal  Article:  ” Metrics  for  Evaluating  Feature-Based  Mapping  Performance ”,  Pablo  Barrios, 
Martin  Adams,  Keith  Leung,  Felipe  Inostroza,  Ghayur  Naqvi,  Marcos  E.  Orchard,  IEEE  Transac¬ 
tions  on  Robotics,  Vol.  33,  Issue  1,  pages  198  to  213,  Feb.  2017. 
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•  Conference  Article:  ”  Estimating  Detection  Statistics  within  a  Bayes- Closed  Multi- Object  Filter”, 
Javier  Correa,  Martin  Adams,  19th  International  Conference  on  Information  Fusion  (FUSION), 
July,  2016. 

•  SSA2016  Workshop  on  Space  Situational  Awareness,  La  Serena,  Chile,  April  2016.  Presentation 
title:  ” Joint  Object  Detection  and  Tracking  Filter  for  Chilbolton  Advanced  Meteorological  Radar 
Data  Processing”. 

3  Datasets 

In  this  section,  a  brief  description  of  the  data  sets  used  is  provided. 

3.1  CAMRa  Radar  Data  sets 

Two  types  of  data  sets  were  obtained  from  CAMRa:  raw  and  post-processed.  For  each  data  acquisition 
step,  the  whole  A-scope  of  a  radar,  namely  the  function  of  returned  signal  intensity  over  the  range,  is 
recorded  and  stored  in  a  raw  (netCDF)  file,  which  is  typical  for  an  atmospheric  radar.  Each  A-scope 
contains  27,000  75-m  bins  for  returned  co-polar  and  cross-polar  signal  strength.  The  size  of  a  raw  data 
file  averages  at  2  Gb  per  observation.  For  size  considerations,  the  raw  data  is  post-processed,  leaving 
one  detection  per  time  step,  associated  with  the  strongest  returned  signal  strength.  For  a  more  detailed 
explanation  of  the  CAMRa  the  viewer  is  referred  to  [5] . 

3.2  The  OmniSSA  Visual  Data  set 

This  OmniSSA  data  sets  consist  of  images  taken  from  three  identical  SR  cameras  positioned  equidistant ly. 
The  camera  resolution  is  3326  x  2054  pixels,  the  field  of  view  is  66°  x  82°.  The  three  cameras  take 
the  pictures  simultaneously,  and  their  calibration  parameters  using  the  Brown  distortion  model  [6]  are 
included.  The  data  set  is  composed  of  observations  made  during  one  night  in  Atlanta  (33.777468°  N, 
84.398969°  W)  and  it  contains  27  images  per  camera  in  the  Flexible  Image  Transport  System  (FITS) 
format.  The  images  were  taken  at  different  exposition  times,  the  first  1  sec.,  thee  of  them  15  sec.,  and  the 
rest  30  sec.  The  database  was  provided  by  Dr.  Marcus  Holzinger  from  Georgia  Institute  of  Technology. 
Further  details  about  the  data  set  can  be  found  in  [7]. 

Figure  1  shows  the  sampling  times  of  the  data  set.  The  vertical  axis  is  the  label  of  the  image. 
The  horizontal  axis  is  the  time  the  photographs  were  taken,  while  the  length  of  the  line  represents  its 
exposition  time.  Figure  lb  shows  a  subsampling  between  images  7  to  11.  Figures  show  that  the  sample 
time  is  irregular,  and  at  the  best  case  the  same  object  (satellite)  appears  in  two  following  images.  For 
this  reason,  this  data  set  cannot  be  used  for  tracking,  but  other  tasks  can  be  done,  like  object  detection 
or  coordinate  mapping  (conversion  from  image  pixels  to  celestial  coordinates  and  vice  versa). 


4  Radar-based  Tracking 

In  this  section,  methods  two  methods  for  processing  the  data  generated  by  the  radar-based  tracker  are 
presented.  For  a  detailed  explanation,  please  refer  to  the  attached  article. 

4.1  Constant  False  Alarm  Rate  (CFAR)  Processor 

The  original  detection  extraction  algorithm  in  CAMRa  measurement  processing  was  using  the  range  bin 
corresponding  to  the  strongest  signal  strength.  This  method  could  fail  to  detect  a  valid  target  in  the 
presence  of  occludents,  transient  objects  or  high  intensity  noise.  A  more  robust  method  of  line  of  sight 
peak  detection  is  the  CFAR  family  of  algorithms  [1].  The  family  of  CFAR  algorithms  are  based  on  local 
statistics  estimation  to  estimate  peaks  within  the  received  radar  signal.  In  particular,  the  CA-CFAR 
algorithm,  diagrammatically  shown  in  Figure  2,  is  evaluated  in  this  report.  The  CA-CFAR  algorithm 
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Figure  1:  Sample  times  of  the  data  set 


can  be  described  as  a  sliding  window  operator  around  the  current  cell  that  evaluates  its  neighborhood 
to  determine  the  local  threshold.  If  the  current  cell  being  evaluated  is  above  the  local  threshold,  it  is 
reported  as  a  peak.  Since  this  method  works  on  local  basis,  it  can  report  more  than  one  peak  per  radar 
scan. 


Detection 


i  i  CUT  i  i 

v - V - ' 

Inner  window 

v - v - y 

Outer  window 

Figure  2:  An  example  of  theCA-CFAR  operation.  The  neighborhood  of  the  current  cell  under  test  (CUT) 
is  evaluated  to  estimate  an  appropriate  detection  threshold.  The  black  square  designates  the  CUT  and 
the  corresponding  guard  cells  are  shown  in  Grey. 


4.2  Gaussian  Mixture  Joint  Target  Detection  and  Tracking  Filter 

For  the  original  analytical  explanation  of  the  JoTT  filter  the  reader  is  referred  to  [21,  12,  16]. 

The  JoTT  filter  is  a  RFS-based  filter  designed  to  estimate  a  single  track  when  multiple  measurements 
are  received.  To  do  so,  the  presence  of  an  object  and  its  state  is  modeled  as  a  Bernoulli  RFS: 


X  -  Bernoulli  RFS  (q,  S )  =  /Bernoulli  (X) 


fl-q ,  if  X  =  0 

<  qs(x),  if  X  =  {x} 
0,  otherwise, 


(i) 


with  X  being  a  random  set  that  can  one  zero  element  with  probability  1  —  q  and  exactly  one  element 
with  probability  q. 

Analogous  to  standard  vector-based  filtering,  a  set-transition  model  is  required.  A  standard,  or 
“natural”,  transition  model  for  an  RFS  state  is  shown  in  Figure  3  and  defined  as: 

•  If  no  target  is  present,  one  can  be  born  with  probability  Pg,  and  distributed  according  to  /#(•). 
Conversely,  if  no  target  is  present,  there  will  be  no  target  with  probability  1  —  Pb-  These  cases  are 
shown  in  Figure  3a. 
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•  If  a  target  is  present  in  the  environment,  the  target  survives  with  probability  Ps(x)  and  its  state 
evolves  in  time  according  to  fk\k-i('\')-  If  a  target  is  present,  it  will  cease  to  exists  with  probability 
1  —  Ps-  These  cases  are  shown  in  Figure  3b. 

Using  this  description,  the  pdf  for  the  transition  RFS  P(Xk \Xk-i)  can  be  expressed  as  follows: 


P{Xk\Xk_1) 


Bernoulli  RFS  (Ps(xfe_i),/fe|fe_1(-|xfc_i)),  if  Xfe_i  =  {xfe_i} 
Bernoulli  RFS  (PB,  fB),  if  Vt-i  =  0, 


(2) 


where  Ps(x)  is  the  probability  of  survival  of  a  target  located  at  x,  /fc|fc-i(xfc|xfc_i)  is  the  prediction 
model  for  a  target  located  at  x^-i,  Pb  is  the  probability  of  the  target  appearing  in  the  environment  and 
/b(x)  is  the  spatial  distribution  of  the  newly  born  target. 

Now,  assuming  that  P(Xk-i\Zi:k~i)  is  Bernoulli  distributed  -  i.e.  Bernoulli  RFS  (qk- 1,  s^-i),  it  is 
possible  to  show  that  the  predicted  RFS,  P(Xk\Zi:k~i),  is  also  a  Bernoulli  RFS  [12,  Section  14.7.3]  with 
parameters: 


Qk\k-i 
$k\k—l  (Xfc) 


(1  —  Qk-l)PB  +  Qk- 1  Ps 

(1  -gfe-l)^fe(Xfc)  Pqk-lPsi^k)  f  fk\k-l(xkjxk-l)sk-1(xk-1)dxk-1 

Qkjk-1 


(3) 

(4) 


The  predicted  RFS,  can  be  corrected  with  the  received  set-based  measurements  using  Bayes  theorem 
to  obtain  the  required  estimated  posterior  RFS.  To  achieve  this,  first  a  model  of  the  measurements  is 
required.  Similar  to  the  standard  multi-target  model,  the  standard  multi-target  detection  model,  shown 
in  Figure  4,  is  defined  as: 


1.  Received  measurements,  shown  as  squares  in  Figure  4,  are  either  clutter  or  a  single  measurement 
resulting  from  the  target,  shown  as  circles  in  Figure  4. 

2.  If  a  target  is  present  in  the  environment,  it  generates  a  measurement  with  probability  Pd(x)  which 
is  distributed  according  to  fz( z|x).  This  case  is  shown  in  Figure  4a. 

3.  Clutter  measurements  are  statistically  independent  of  the  target. 


This  model  can  be  interpreted  as  the  union  of  two  sets, 


Z  =  O  U  Z\X,  (5) 

namely  the  set  of  clutter  measurements  0  and  the  set  with  the  valid  measurement  generated  by  the 
target  Z\X.  A  standard  model  for  the  clutter  RFS  0  is  the  Poisson  RFS  with  parameters  A©  for  the 
clutter  rate  and  k(z)  for  the  spatial  distribution  of  the  clutter  measurements,  as  shown  in  Figure  4b.  In  a 


1  -Pb 


1  -Ps 


k  —  1 


k-1 


•  ~  /s(Xfc) 


•  ~  /fc|fc-l(Xt|Xfc_i) 


Pb 


(a)  Transition  model  for  Xk-±  =0 

(b)  transition  model  lor  Xk—i  —  {x^—i} 
Figure  3:  Graphical  representation  of  the  standard  RFS  transition  model. 
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X 


1  -Pd 


X 


fz{z  |x) 


z 

K©(zi) 

K@(Zn) 


Poisson(A@) 


Pd 


(b)  Observation  model  for  clutter  measure- 
(a)  Observation  model  for  the  target  X  ments 


Figure  4:  Graphical  representation  of  the  standard  RFS  detection  model. 


Poisson  RFS  the  number  of  elements  of  the  set  is  Poisson  distributed,  while  the  elements  are  identically 
and  independently  distributed  according  to  the  spatial  distribution. 

Following  the  description  of  the  standard  multi-target  measurement  model,  the  model  for  the  mea¬ 
surement  resulting  from  the  target  can  be  written  as: 


P(Z\X) 


Bernoulli  RFS  (Pd  (x),  fz(z  |  x)),  if  X  =  {x} 

0,  otherwise 


(6) 


where  fz( z|x)  is  the  likelihood  of  observing  z  from  target  state  x  and  Pd(x)  is  the  probability  of  detecting 
a  target  located  at  x. 

It  is  possible  to  show  that,  using  the  standard  multi-target  likelihood,  if  the  prior  distribution  is 
distributed  as  Bernoulli  RFS  (qk\k-h  fk\k-i)  then  the  posterior  distribution  is  also  a  Bernoulli  RFS  [12, 
Section  14.7.4]  now  with  parameters: 


/  (1  -  PD(x))sft|fe_1(x)  +  EzeZfe 


Qk  —  Qk\k-1 - 7 - 

1  -  qk\k-§  +  <lk\k-i  [f  (1  -  fb(x))sfe|fe_1(x)  +  Ez ez 


Aqk(z) 

Pp(x)/Z(z|x)sfe|fe_1(x)  ^ 


A©/^(z) 


Sfc(x)  =  —  [  (1  —  PD(x))sfe|fc_1(x)  +  y] 
Qk  \ 


-Pp(x)E(z|x)sfc|fc_1(x) 
Aqk(z) 


(7) 

(8) 


An  analytic  solution  for  this  filter  is  obtained  when  the  underlying  spatial  distribution  s(x)  is  assumed 
to  be  a  Gaussian  mixture: 


N 

Skfe)  =  ^,%(x,mW,pW).  (9) 

i= 1 

Similarly,  the  distribution  of  the  locations  where  the  target  can  appear,  is  also  modeled  with  a  Gaussian 
mixture: 


N 

/b(x)  =  X.wbX(^^>,P(s>),  (10) 

i= 1 

where  Wg\  and  are  the  weights,  means  and  covariances  of  the  components  of  the  birth  process. 
The  the  probabilities  of  survival  and  detection  are  constant  in  this  formulation  of  the  filter. 

With  these  assumptions  it  is  possible  to  show  that  after  the  prediction  step,  the  new  parameters  of 

10 


DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 


the  Bernoulli  RFS  are: 


Qk\k-i  =Pb{  1  -  qk)  +  PsQk 

/  \  Pb{  1  ~  qk)  r  /  x  , 
Sfc|fe-iW  = — — - /sfcw+ 


Nk-! 


PsQk 


Qk\k-1  i=1 


4-i^(x;mSfe-i;pfc|Li)- 


(ii) 


(12) 


The  distribution  of  the  components  of  the  mixture  are  also  preserved  as  Gaussians,  with  mean  m/c|/c_1 
and  covariance  P/C|/C_i  corresponding  to  the  regular  Kalman  filter’s  or  Extended  Kalman  Filter  (EKF)’s 
prediction  equations: 


mk\k-i  =Fjfe_im^| 


1/c- 

>w 


(13) 

P  fe|fc-i  =F,_iP-iFLi  +  Q/c— i  •  (14) 

During  the  update  step  of  the  JoTT  filter,  the  corrected  parameters  for  the  Bernoulli  RFS  are: 

1  —  Afc 


qk 


1  Qk \k  —  1  A/e 

1  ~Pd 


Qk\k—1 


^  1  —  Afe 

r 

lAfc  =Pd 


iVfc_i  (j)  (i),  , 

5/c|/c-i(x)  +  :  _  A  \  A/ (x,  mfc  ,  F*.  ) 


zCZ  2=1 


\k(z) 


Nfc-1 

,  \"  \'  wfc  -^fc  tV 

A-  A?  Ak(z) 

ZG^  2=1  '  ' 


(15) 

(16) 

(17) 


where: 


>(z)  tW( z;Hftm«_1;H,P«_1Hr+Rfc) 

(18) 

Kfc  =Pgc_1Hj+Rfe 

(19) 

mfc  =mfc|fe_1  +  (z  - 

(20) 

P/e  — (I  -  KkHk)I>k\k-l- 

(21) 

4.3  Robust  Versions  of  the  JoTT  Filter 

The  standard  JoTT  filter,  requires  the  parameters  for  the  detection  statistics  to  be  specified  a  priori, 
namely  the  probability  of  detecting  the  target  and  the  rate  at  which  clutter  is  detected  and  how  it  is 
distributed  spatially.  These  parameters  are  non-trivial  to  describe  in  terms  of  a  theoretical  analysis  of 
the  sensor  or  to  estimate  in  an  off-line  process,  since  the  detection  statistics  are  likely  to  vary  in  time. 
To  overcome  this  issue,  two  extensions  to  the  JoTT  filter  were  developed  which  are  capable  of  jointly 
estimating  the  presence  or  absence  of,  and  the  distribution  of  a  target  and  the  probability  of  detection 
and  clutter  rate. 

4.3.1  JoTT  filter  with  unknown  Pd 

In  the  formulation  of  the  JoTT  filter,  it  was  assumed  that  the  probability  of  detection  was  state  dependent. 
This  suggests  that  we  could  extend  our  state  space  to  directly  include  the  probability  of  detection  as  one 
its  values, 

x=[x,Pd]-  (22) 

Given  that  a  model  for  the  distribution  of  Pd  is  provided,  this  simple  extension  allows  the  filter  to  include 
the  estimation  of  Pd  at  a  very  low  cost. 
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To  account  for  a  possibly  unknown  and  time  varying  Pd,  a  /3GM  was  used  in  [13]  as  a  variation  of 
a  Gaussian  Mixture  for  the  Po-Probability  Hypothesis  Density  (PHD)  and  Pp-Cardinalized  Probability 
Hypothesis  Density  (CPHD)  filters. 

The  appearance  or  disappearance  of  the  target  in  the  environment  can  be  seen  as  a  Bernoulli  exper¬ 
iment,  and  a  natural  model  for  the  parameter  of  a  Bernoulli  distribution,  in  this  case  the  probability  of 
detection,  is  a  Beta  distribution  [4,  Section  8.2].  As  such,  and  similar  to  the  standard  JoTT  filter,  an 
analytical  solution  is  obtained  assuming  that  the  distribution  of  the  new  state  shown  in  Equation  22  is 
now  a  Beta-Gaussian  mixture. 

The  prediction  step  for  the  probability  of  target  existence  and  the  spatial  distribution  is  similar  to 
the  standard  JoTT  filter.  They  differ  in  that  the  spatial  birth  distribution  has  to  have  Beta  components 
along  with  the  Gaussian  component: 

N 

/b(x)  =  JV (x,  mg  ,Pg  )  Beta(P/5,  &g*).  (23) 

i=l 


The  prediction  model  for  the  Beta  assumes  that  the  probability  of  detection  follows  a  random  walk. 
In  practice  this  is  implemented  by,  at  the  prediction  step,  maintaining  the  expected  value  and  increasing 
the  variance  by  a  small  amount.  This  avoids  the  problem  of  the  Beta  components  converging  to  a  highly 
concentrated  point,  and  allows  the  incorporation  of  temporal  variability.  This  is  achieved  by  multiplying 
both  parameters  of  the  Beta  distribution  by  a  constant  value,  resulting  in  the  maintenance  of  the  expected 
value,  but  an  increase  in  the  variance  by  a  value  e.  The  multiplicative  constant  can  be  estimated  from: 


a 


2  (i) 

B,k\k—1 


=  < 7 


2  (i) 

B,k-1 


+  e. 


(24) 


This  results  in  the  following  parameters  for  the  predicted  Beta  components 


Q{i) 


ak\k—l 
"a- |  A-  :i 


i(i)  6(i)  —  e 


i. 


(25) 

(26) 
(27) 


Once  the  predicted  distribution  is  obtained,  the  ydGM-JoTT  update  resembles  the  JoTT  update  with 
the  slight  difference  that  each  instance  of  the  Pd  term,  which  appears  in  the  original  JoTT  filter,  is 

a(i)  _  b{i)  _ 

replaced  by  (i)  klk  ^ — .  Conversely,  each  instance  of  the  term  (1  —  Pd)  is  replaced  by  (i)  fc|fc  m — • 

ak\k-i+bk\k-i  ak\k-i+bk\k-i 

This  results  in  Equation  16  being  replaced  by 


Nt 


k\k-l 


$k  ([x,  PD}) 


1  —  At 


wk\k~i 


b{i) 

uk\k~i 


i=  1 


k\k—l  ^  uk\k—l 


V  (x;  Pfc|i-i)  Beta(Po,  a$,fc>  &$*) 


.  --k\k—  1  n(i)  T  (i)  (  \ 

I  1  ^1^-1  Lk  \z)  Kr(„.  mW.pW^RrtnfP  n®  h®  i  focA 

+  l  —  /\  Z_^  (i)  Ai)  Xk(z)  ^  X’  m/c  aD,k'  °D,kJ’ 


„  .  ,  /)-  _I_  hyi}  Att(z) 

zez  i=l  ak\k-l  '  °k\k-l  V  ’ 

with  the  new  Beta  components  (assuming  i  —  1 . . .  Nk\k_i  and  j  =  1 . . .  \Z\) 


( i )  _  (i) 

aM,k  —  ak\k-l’ 

(Nk\k-i+ij)  _  (i)  I  -I  ANk\k-i+ij)  _  7  (i) 

aD,k  ~ak\k-l^Li  UD,k  ~  Uk\k-V 


b{i) 

uM,k 


fob)  4-  1 
°k\k-l  +  1 
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4.3.2  JoTT  filter  with  unknown  detection  statistics 

The  formulation  presented  in  the  previous  section  has  its  limitations.  For  example,  it  is  not  possible  to 
extend  the  state  space  to  account  for  clutter,  since  in  the  derivation  of  the  JoTT  filter,  clutter  and  targets 
are  assumed  independent  from  one  another,  and  thus  the  clutter  parameter  would  be  unobserved. 

To  estimate  the  clutter  process,  Mahler  [13,  Section  18.3]  used  the  idea  of  clutter  generators.  Clutter 
generators  are  analogous  to  targets,  but  generate  clutter  measurements  only.  By  applying  this  model  to 
the  JoTT  filter,  it  would  convert  a  single-target  filter,  into  a  multi-target  filter,  as  it  would  be  required 
to  keep  track  of  the  object  of  interest  as  well  as  the  possible  multiple  clutter  generators.  A  different 
approach  is  used  in  this  work,  where  a  standard  Poisson  model  is  assumed  for  the  clutter  process,  but 
with  an  unknown,  and  possibly  time- varying,  clutter  rate  parameter.  The  problem  then  becomes  how  to 
model  and  estimate  the  clutter  rate  parameter. 

To  derive  a  JoTT  filter  robust  to  both,  unknown  probability  of  detection  Pd  and  clutter  rate  A©,  a 
full  probabilistic  re-formulation  is  required.  Instead  of  extending  the  state  space,  both  these  statistics  are 
incorporated  into  the  estimation  problem.  Thus,  the  values  which  must  be  estimated  are  the  Bernoulli 
set  X  corresponding  to  the  estimated  target,  the  probability  of  detection  Pd  and  the  rate  at  which  clutter 
appears  in  the  sensor  data  A©,  resulting  in  a  joint  probability  distribution 


P(Xk,PDk,  A©JZi:fc), 


(29) 


with  Xk  being  the  set  representation  of  the  target  (empty  or  single  element  set),  is  the  probability 
of  detection  and  A©fc  is  the  clutter  rate  at  time  step  k. 

This  model  has  a  secondary  advantage.  Since  the  value  for  Pd  does  not  depend  on  the  estimated 
state,  it  is  now  theoretically  possible  to  obtain  and  estimate  Pd  even  when  no  target  is  present  in  the 
environment  (e.g.  q  m  0).  This  helps  in  the  (re-) initialization  of  targets  when  there  has  been  a  long 
period  of  only  clutter  measurements. 

It  is  assumed  that  the  joint  probability  distribution  in  (29)  has  the  following  form: 


P(Xk,PDk,\Qk\Zl:k) 


'  Vm  «;(i)  <7(i)  R(i) 

2^=1  wQ,kuQ,knQ,k 

vra  ■>„(*)  ru>  rW  „(») 

2^i=l  wX,k^rX,k1-,X,khX,k 


,  if  X  =  0 
,  if  X  =  {x}  , 
,  otherwise, 


(30) 


with  Gfjk  =  Gfjk(\ ©),£>□£  =  Ba\(PD)  and  s^k  =  sn^(x)  being  the  distributions  for  the  clutter 
rate,  probability  of  detection  and  the  target’s  spatial  distribution  respectively,  and  w^k  the  weights 

such  that  w0k  wx\  =  1-  The  square  symbol  can  be  either  □  =  0,  belonging  to  the 

mixture  representing  target  non-existence,  or  □  =  X,  indicating  target  existence  in  Equation  (30).  The 
arguments  of  the  functions  GjjJk,  B ^  k  and  Sk  will  be  omitted  when  possible  for  a  clearer  exposition.  For 
a  fully  robust  filter,  with  unknown  detection  statistics,  the  complexity  of  the  formulation  of  the  filter 
increases  considerably.  Intuitively,  the  probability  of  a  target  being  present  in  the  environment  is  now 
q  =  i  wx\>  while  the  probability  of  no  target  being  present  is  1  —  q  =  Xq=i 

Assuming  that  the  number  of  clutter  measurements  is  Poisson  distributed,  a  natural  distribution  to 
model  the  rate  of  clutter  is  the  Gamma  distribution  [4,  Section  10.1].  As  such,  a  Gamma  distribution  is 
adopted  to  model  the  clutter  rate, 


Gc u(V>)  =  Gamma(Ae,a^y /3^fe) 


A 


w 


AO 

*U,k 


r(«gfc) 


*<0 
\  a, k 

“A© 


-l 


e-O0 


(31) 


(£)  f£\ 

where  c^k  and  P^k  are  the  parameters  of  the  Gamma  distribution,  and  a  Beta  distribution  to  model 
the  detection  probability, 


AO 


-l , 


AO 


Bn,k(pD)  =Beta(PD,a£> 


W  hW  \  _ 


D 


(1  -Pd)°^ 


-l 


P(aa,k 


bw  ) 


(32) 
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with  a^jk  and  b^k  being  the  parameters  of  the  Beta  distribution. 

At  the  prediction  step,  assuming  that  the  distribution  at  time  k  —  1  is  of  the  form  of  Equation  (30), 
it  is  possible  to  show  that  the  predictive  distribution  has  the  same  form: 


( V  R0) 

p/Y  p  \  |7  \  _  J  ^j=l  W0,fc|fc-l'-T0,fc|fc-|-D0,fc|fc-l 

-L\^k\k— 1 5  -PD/s|fe-i  5  ,1  — l)  \  ^  (i)  o(d  (d 

[2^i=  l  wX,k\k-lLrX,k\k-l^X,k\k-lSX,k\k-l 

,  if  —  0 

,  if  X  =  {x}, 

(33) 

where  for  the  first  m  elements  of  the  mixture  (i  =  1 . . .  m) 

w(X,k\k-l  =  ^0,1-1^’  W^k\k-1  =  ^0,1-1  (1  _ 

(34) 

Ad  _ /Ad  Ad  _ Ad 

°X,/c|/c-l  —  ^0,^-1’ °'0,/c|/c-l  —  'Jr0,fc  — 1 5 

(35) 

AO  _  o(d  AO  _  o(d 

(36) 

(d  f 

5X,/c|/c-l  —  4b, 

(37) 

and  for  the  last  n  elements 

dbifc-i  =  Afcife-i  =  -  **?), 

(38) 

/Ad  _  AO  A*+m)  _  AO 

Lt0,/c|/c-1  — 

(39) 

p(i+m)  _  p(d  p(*+m)  _  p(d 

_  n0,k\k-i  —  nx,k-n 

(40) 

42-1  =  ~  /  -ps(xfc-i)/fc|fc-i(xfc|fc-i|xfc_i)s^/s_;|(xfe_i)cbcfc_ 
Tz  J 

-i5 

(41) 

and  Ti  is  a  normalization  constant  calculated  as  follows 


Ti  =  JJ  Ps(xfe_i)/fc|fc_i(xfe|fe_1|xfc_i)s^fc(xfc_i)dxfc_idxfe.  (42) 

Note  that  the  parameters  related  to  the  clutter  and  probability  of  detection,  under  this  model,  should  re¬ 
main  constant.  Similar  to  the  /3GM-JoTT  filter,  to  avoid  overconfidence  and  allow  a  small  time  variation, 
a  random  walk  is  used  to  model  the  evolution  of  both  the  clutter  rate  and  probability  of  detection. 

For  a  clearer  exposition,  we  will  assume  a  constant  probability  of  survival,  Ps(’)  =  -Ps,  which  then 
results  in 


^  =  Ps  (43) 

SX,k(Xk\k-l)  =  J  fk\k^i(^k\k-l\^k~rl)Sxtk-l(Xk-l)dxk-l.  (44) 

To  account  for  the  observation  at  the  current  time  step,  if  the  prior  distribution  is  assumed  to  be  of 
the  form  shown  in  Equation  (30),  by  using  Bayes  theorem  with  the  standard  multi-target  likelihood  it 
can  also  be  shown  that  the  corrected  (or  posterior)  probability  distribution  has  the  same  form 


P(Xk,PDk,\ek\Z1:k) 


E 

E 


n/j)  r<{i)  pU) 

■j= 1  w0,ku0,kn<D,k 


m(l  +  |Z|)„„(i) 
i=l 


Ad  Ad  „(d 


wy~'  } 

w  X ,  k  ^  X ,  k  D  X ,  k b  X ,  k 


,  if  Xk  =  0 
,  if  Xk  =  {x},  ’ 


where,  for  the  case  of  Xk  =  0,  the  parameters  of  the  mixture  are 


(45) 


"  0.A-  —  "  0.A-  /,•••  ! 


r  [a 


,U) 

\,k\k-l 


:\)0, 


,  a(j) 

(j)  a0,fc|fc-l 

k\k  —  l 


(j) 

fa(j)  )  f/3(j)  +1)a0-klk 

{a0,klk-l  J  l^0,fc|fc-l  +  1J 


a 


U) 


k  ~  a®%k-l  +  1^1  ~  PQ,k\k-l  + 


U) 


),k\k- 


(46) 

(47) 
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and  the  parameters  of  the  Beta  component  B^k  are  the  same  as  the  B^k_ r  Now,  for  the  mixture  in  the 
case  of  Xk  =  {x},  we  have  m  components,  related  to  miss-detecting  the  target,  with  parameters 


WX,k  —  WX,k\k-l  (i) 


b(i) 

uX,k\k-l 


r  (a 


,w 

'X,k\k  —  1 


(i) 

( i )  aX,k\k-l 

X,k\k-1 


a 


x]k  ~  a(X,k\k-l  +  \Z\>0X,k  ~  Px]k\k-1  +  ^ 


o(0  -0w  fc(0  -1,(0  +1 

aJC,/c  —  aX,k\k-V  UX,k  —  uX,k\k-l  +  x> 

and  m\Z\  components,  related  to  detecting  the  target,  with  parameters  (j  =  1 . .  .to, £  =  1 . . .  |Z|) 


(48) 

(49) 
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4.4  State  and  uncertainty  propagation. 

Each  Gaussian  component’s  mean  and  covariance  E/c|/c_1  in  a  mixture  is  managed  using  the 

corresponding  Unscented  Kalman  Filter  (UKF)  equations  [10].  The  propagation  of  the  target  state 
and  sigma  points  is  performed  through  the  numerical  integration  of  the  state  vector  according  to  the 
model  where  only  gravitational  interaction  between  two  bodies  is  taken  into  account.  In  this  model,  the 
gravitational  acceleration  ac  is  represented  as: 

aG  =  ~W ’  (56) 

where  /jl  is  the  Earth’s  gravitational  constant.  The  derivative  of  the  state  vector  x  used  in  the  integration 
is  then: 


dt 


X4  £5  Xq 


~l-LXl 


—  fj-x  2 


—yxi 


1  T 


(V xl+xl+xl)3  (yTf+Tf+Tf)3  xl+xl+xl)3  J 


(57) 


5  Vision-based  Tracking 

Initial  work  has  also  been  carried  out  using  images  for  detecting  and  tracking  space  debris.  This  section 
describes  the  work  carried  out  to  analyze  the  OmniSSA  data  set  in  order  to  detect  possible  objects  of 
interest.  The  OmniSSA  data  set  does  not  contain  data  for  performing  tracking,  but  the  images  can  be 
used  for  pre-processing,  such  as  reducing  background  noise,  removing  stars  and  other  objects  from  the 
image  to  be  analysed,  object  detection  and  obtaining  a  mapping  between  celestial  coordinates  and  the 
image  pixels. 
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5.1  Object  detection 

SExtractor  [3]  is  a  widely  used  software  for  source  extraction  from  astronomical  images.  It  detects  objects 
such  as  stars,  satellites,  galaxies  from  FITS  images.  Then  it  computes  photometry1  from  the  detected 
objects  and  creates  a  catalog  with  information  including  position,  dimension,  shape  and  their  statistics. 

Figure  5  shows  a  diagram  of  the  main  stages  that  SExtractor  performs  for  building  a  catalog  of 
objects  from  an  astronomical  image.  The  process  starts  by  estimating  the  background  produced  by  noise 
presented  in  the  image.  The  image  is  divided  in  blocks  and  a  local  background  noise  estimate  is  obtained 
from  each  block  using  an  iterative  method.  For  each  iteration,  the  median  ( block )  and  standard  deviation 
((j block)  of  the  pixel  intensities  are  computed.  Pixel  intensities,  which  differ  significantly  from  the  median 
value,  i.e.  those  that  fulfill  | I(x,y)  —  hiockl  >  3cr  are  rejected,  and  a  new  iteration  starts  with  the 
remaining  pixels.  The  variables  x  and  y  indicate  the  pixel  position  in  the  image.  Iteration  stops  when  no 
pixels  are  rejected.  The  mean  of  the  remaining  pixels  is  considered  as  the  block  background  value,  unless 
the  (Jbiock  value  changes  more  than  20%,  in  which  case  the  field  is  considered  crowded  and  the  mode  value 
is  taken  instead.  After  computing  the  background  for  all  blocks,  a  median  filter  is  applied  over  the  entire 
image  in  order  to  eliminate  possible  local  extreme  intensities  produced  by  bright  stars.  As  the  image  was 
divided  in  blocks,  the  temporal  background  image  is  low-resolution,  thus  in  order  to  obtain  the  full-size 
background  image  /#(£,?/),  a  bicubic-spline  interpolation  is  applied  to  the  low-resolution  image. 

The  background-subtracted  image  Ip(x,y)  is  the  difference  between  the  input  image  and  the  back¬ 
ground  image.  There  is  an  option  of  using  a  filter  in  order  to  smooth  Ip(x,y).  The  image  Ip^x^y)  is 
convolved  with  the  filter,  helping  to  improve  the  detection  of  faint  and  extended  objects.  SExtractor 
provides  four  types  of  filter:  Gaussian,  Mexican- hat,  top-hat  and  block,  all  of  them  normalized  and  at 
various  sizes,  although  other  filters  can  be  applied  [2,  9].  Gaussian  (Fg(x,t/)),  top-hat  ( Fth(x,y ))  and 
block  (Fb(x,y))  filters  compute  a  weighted  average  around  the  pixel,  useful  for  emphasizing  object  pixels 
while  reducing  the  surrounding  noise.  Mexican- hat  (FMh{x,y))  is  a  passband- filter,  which  smooths  the 
image,  but  also  emphasizes  edges.  It  is  useful  in  crowded  star  fields.  The  filter  functions  are  shown  in 
Equations  58.  For  Fth(x,y),  r  represents  the  radius  of  the  ”  top-hat”.  a2  represents  the  variance  of  the 
Gaussian  part  of  Fc(x,y)  and  Fmh{x^v)- 


Fc(x,y) 

Fth(x,y) 

Fb(x,y) 

FMh{x,y) 


1 

27 rcr2 


exp 


(  X2  +y2\ 

V  2,2  j- 


1,  if  x2  +  y2  ^  r2, 
0,  otherwise. 


x2  +y2\ 
2<t2  ) 


exp 


(58) 


A  threshold  (th)  is  applied  to  the  image  and  those  pixels  for  which  Ip(x,y)  >  t h  are  considered  to  be 
part  of  an  object.  A  deblending  algorithm  is  used  to  split  blobs  formed  by  more  than  one  object.  For 
details  of  the  deblending  algorithm  refer  to  the  manual  [2].  The  segmented  image  is  denoted  Is{x,y). 
Finally  the  catalog  of  astronomical  objects  is  built  by  applying  photometry  to  the  segmented  objects. 

Other  useful  images  can  be  computed  by  SExtractor,  and  used  for  posible  future  debris  detection 
algorithms,  such  as  a  background-subtracted  image  in  which  stars  are  removed  (/r(x,?/)).  For  example 
debris  can  be  identified  as  straight  lines  with  pixel  intensity  values  higher  than  the  background,  but  lower 
than  the  stars.  This  occurs  because  stars  produce  their  own  light,  whereas  orbiting  elements  do  not.  In 
order  to  get  /#(%,?/)  from  SExtractor,  a  suitable  threshold,  higher  than  the  intensity  of  the  debris  but 
lower  than  stars,  must  be  selected. 

The  next  Section  explains  how  star  detection  is  used  to  obtain  the  transformation  between  image 
coordinates  and  celestial  coordinates. 

1  Cited  from  [17].  Photometry  is  the  science  of  measuring  the  flux  received  from  celestial  objects.  It  usually  refers  to 
measurements  of  flux  over  broad  wavelength  bands  of  radiation. 
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Figure  5:  Flow  diagram  of  the  most  important  components  of  SExtractor  showing  how  it  builds  a  catalog 
of  objects  from  an  astronomical  image. 

5.2  Coordinate  Extraction  from  the  Images 

In  order  to  map  pixel  coordinates  from  the  images  to  celestial  coordinates  the  SExtractor  [3]  and  As- 
trometry.net  [11]  software  are  used.  The  first  step  is  to  undistort  the  original  images  using  the  provided 
calibration  parameters,  then  stars  and  other  objects  are  detected  with  SExtractor,  and  finally  the  detected 
stars  are  compared,  matched  and  identified  with  a  catalog  of  stars.  The  Tycho-2  catalog  [8],  provided  by 
Astrometry.net,  is  chosen  for  this  task  because  it  works  better  for  wide-angle  images.  The  result  of  the 
mapping  is  a  polynomial  that  maps  pixels  to  right  ascension  and  declination  from  the  equatorial  celestial 
coordinates,  following  the  World  Coordinate  Systems  (WCS)  standard.  The  equatorial  coordinates  given 
by  Astrometry.net  are  topocentric,  and  the  camera  position,  latitude  and  longitude,  is  needed  to  obtain 
the  geocentric  coordinates. 

Debris  identification  can  be  carried  out  by  extracting  known  objects  from  TLE  catalogs.  Debris 
locations  are  predicted  and  then  projected  into  the  image  space  using  the  previously  obtained  mapping. 
Details  and  some  examples  can  be  found  in  the  Results  section. 


6  Results 

In  this  section  results  of  the  developed  robustified  single  target  multiple  measurements  filters  and  of  the 
vision-based  detection  algorithms  are  presented.  The  JoTT  filter  variants  are  evaluated  on  the  radar- 
based  CAMRa  data  set,  while  the  visual  detection  are  evaluated  on  the  OmniSSA  data  set. 

6.1  Radar-Based  Results 

To  evaluate  the  performance  of  the  proposed  extensions  to  the  standard  JoTT  filter,  first,  simulations  were 
carried  out  to  then  use  the  real-world  CAMRa  data  set.  In  the  simulated  environment,  full  knowledge  of 
the  ground  truth  is  available,  thus  allowing  to  accurately  measure  the  estimation  errors  of  the  different 
evaluated  approaches. 

6.1.1  Simulated  Radar  Measurement 

To  show  that  the  proposed  filters  are  robust  to  unknown  and  time  varying  detection  statistics,  a  pass  of  an 
orbiting  object  is  simulated.  The  simulation  was  created  based  on  the  Cosmos- 1544  TLE  and  Simplified 
General  Perturbations  -  4  (SGP4)  propagation  model  [19].  The  rate  at  which  clutter  measurements 
are  generated  is  arbitrarily  assumed  to  vary  sinusoidally  with  time  with  values  between  4  and  8  clutter 
measurements  per  scan,  while  the  probability  of  detection  is  arbitrarily  assumed  to  follow  a  bell-shaped 
curve,  with  a  minimum  value  of  5%  and  a  maximum  value  of  80%  as  shown  in  Figure  6.  To  process  this 
simulated  data  set,  the  filters  are  not  provided  with  the  real  information  related  to  the  detection  statistics. 
For  benchmark  comparisons,  all  three  filters  were  executed  using  the  same  parameters:  probability  of 
survival  Ps  set  to  99.99%  and  uniform  clutter  spatial  distribution  k(z)  =  25015.77km-3.  For  the  standard 
JoTT  filter,  Pd  was  set  to  80%  and  for  both,  the  JoTT  and  /3GM-JoTT  filters,  A  was  set  to  4  clutter 
measurements  per  time  step.  The  increase  in  variance  for  the  Beta  components  of  the  /iGM-JoTT  and 
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r-JoTT  filters  was  set  to  10  4  per  time  step  and  the  increase  in  variance  for  the  Gamma  component  of 
the  r-JoTT  filter  was  set  to  10-3  per  time  step. 


1.00 
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0  50  100  150  200  250  300 


Seconds 

Figure  6:  The  simulated  data  set.  The  upper  plot  shows  the  ground  truth  value  for  the  probability  of 
detection  versus  time,  modeled  with  a  bell-shaped  curve,  while  the  lower  plot  shows  the  clutter  parameter 
A  which  varies  sinusoidally  with  time.  The  center  plot  shows  the  simulated  measurements  generated  using 
the  corresponding  detection  statistics. 

Qualitative  results  are  shown  in  Figures  7,  8  and  9.  It  can  be  seen  that  in  the  simulated  scenario, 
the  JoTT  and  /3GM-JoTT  have  some  spurious  tracks.  This  is  due  to  the  filters  removing  the  existing 
track,  due  to  the  unmatched  detection  statistics,  and  creating  a  new  track.  This  behavior  is  highlighted 
by  zooming  at  time  steps  0-100,  as  shown  in  Figures  7b  and  8b. It  can  be  seen  that  the  /3GM-JoTT  filter 
suffers  less  from  this  problem  than  the  JoTT  filter.  On  the  other  hand,  the  r-JoTT  does  not  suffer  from 
this  problem  as  shown  in  Figure  9.  It  also  has  to  be  noted  that  the  JoTT  filter  tends  to  have  a  large 
number  of  discontinuities,  as  shown  in  Figure  7b.  This  is  due  to  the  fact  that,  track  continuity,  measured 
by  the  q  value  from  Equation  1,  depends  on  the  value  of  both  the  probability  of  detection  Pd  and  clutter 
rate  A.  Since  the  JoTT  filter  assumes  fixed  values  for  these  parameters,  it  fails  to  provide  a  continuous 
estimate  of  the  target’s  existence  when  parameters  don’t  match  the  real  values.  As  the  /3GM-JoTT  filter 
is  not  dependent  on  the  choice  of  Pd,  it  suffers  less  from  this  problem  and  produces  a  more  continuous 
track  (see  Figure  8b).  The  best  track  is  produced  by  r-JoTT  filter  as  both  Pd  and  A  are  estimated. 

In  Figures  10  and  11,  the  detection  statistic  estimates  of  the  /3 GM-JoTT  and  r-JoTT  filters  are 
shown.  It  can  be  appreciated  that  in  all  cases  the  true  values  are  correctly  estimated,  albeit  with  small 
perturbations.  The  observed  jitters  on  the  estimates  can  be  explained  as  the  prediction  models  for  the 
detection  statistic  parameters  only  increase  their  variance,  and  no  smoothness  constraints  are  assumed. 

Quantitative  results  are  shown  in  Table  1,  where  the  Mean  Absolute  Error  (MAE)  between  the 
estimated  and  the  ground  truth  states  is  shown  in  the  second  column,  while  in  the  third  column,  the 
percentage  of  the  total  simulated  time  that  the  filters  reported  a  target  is  shown.  It  can  be  appreciated 
that  the  JoTT  filter  does  not  perform  well  in  this  simulation.  This  is  due  to  the  inability  of  the  JoTT 
filter  to  differentiate  between  dense  clutter  and  the  real  target  when  the  filter’s  model  of  the  clutter  is 
incorrect.  The  /iGM-JoTT  filter  somehow  overcomes  this  problem.  This  is  clearly  noted  by  the  fact  that 
the  MAE  decreases  to  4.40km.  Clearly  the  r-JoTT  filter  outperforms  other  methods  as  it  is  capable  of 
estimating  both  the  clutter  density  and  the  probability  of  detecting  a  target.  It  has  to  be  noted  that 
the  high  MAE  values  reported  by  the  JoTT  and  /3GM-JoTT  filters  are  mainly  due  to  the  estimated  false 
positives. 


6.1.2  Results  on  the  CAMRa  Data  set 

In  addition  to  qualitative  analysis  of  the  filtering  results  for  the  CAMRa  data  sets,  it  is  possible  to  compare 
the  performance  to  the  corresponding  SGP4  TLE  propagation.  In  order  to  do  so,  a  historic  TLE  related 
to  the  observation  date  was  retrieved  from  Spacetrack  for  the  Cosmos- 1544  satellite.  Superimposed  plots 
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(a)  Estimation  result  using  the  JoTT  filter. 


(b)  Zoom  to  the  estimate  of  the  JoTT  filter  between 
time  steps  0  to  100. 


Figure  7:  Estimated  results  using  the  standard  JoTT  filter. 
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(a)  Estimation  result  using  the  ^GM-JoTT  filter. 


Measurement 

Estimate 


(b)  Zoom  to  the  estimate  of  the  /3GM-JoTT  filter  be¬ 
tween  time  steps  0  to  100. 


Figure  8:  Estimated  results  using  the  /3GM-JoTT  filter. 


Method 

MAE 

%  of  Time  with  Est. 

JoTT 

6.79  Km 

50.66% 

/?GM-JoTT 

4.40  Km 

77.29% 

r-JoTT 

0.08  Km 

99.74% 

Table  1:  Results  for  the  different  single  tracking  approaches.  In  the  first  column,  the  different  methods 
are  shown,  the  second  column  shows  the  MAE  of  the  estimate,  and  in  the  third  column,  the  percentage 
of  the  time  a  valid  estimate  is  reported  by  the  filter  is  shown. 
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(b)  Zoom  to  the  estimate  of  the  r-JoTT  filter  between 
time  steps  0  to  100. 


(a)  Estimation  result  using  the  r-JoTT  filter. 

Figure  9:  Estimated  results  using  the  r-JoTT  filter. 
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Figure  10:  Pd  estimate  using  the  /3GM-JoTT  filter. 
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(a)  A  estimate  using  the  r-JoTT  filter. 


(b)  Pd  estimate  using  the  r-JoTT  filter. 


Figure  11:  Estimates  of  the  detection  statistics  for  the  simulated  scenario  using  the  r-JoTT  filter. 
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(b)  Plots  of  absolute  errors  of  the  estimates  and  measurements  relative  to  SGP4  propagation. 


Figure  12:  Comparison  of  filtering  results  and  SGP4  propagation 


of  measurements,  estimates  and  propagation  as  well  as  estimate  errors  are  shown  in  Figure  12.  The 
absolute  error  for  azimuth  and  elevation  tends  to  be  zero  when  the  target-generated  measurements  are 
present  in  the  scene.  This  is  corroborated  by  the  fact  that  this  particular  TLE  was  used  to  calculate 
radar  look  angels,  i.e.  aforementioned  azimuth  and  elevation.  The  same  kind  of  behavior  cannot  be 
observed  for  the  range  component  as  it  is  known  that  TLE  does  not  provide  precise  position  of  the  object 
of  interest  and  the  father  from  epoch  it  gets,  the  bigger  this  position  error  becomes;  thus,  this  is  the  main 
reason  why  TLE  cannot  be  used  to  generate  reliable  ground  truth. 

Results  of  analyzing  the  multiple- measurement  CAMRa  data  sets  are  shown  in  Figures  13,  14  and  15. 
For  the  single- measurement  filtering  results  please  refer  to  [14]. 

To  test  the  filter  performance  in  real  multiple- measurement  scenario,  two  data  sets  corresponding  to 
two  consecutive  passes  of  the  GPM  satellite,  were  chosen.  Figure  13  shows  the  JoTT  filter  estimates. 
The  parameters  used  to  obtain  these  result  were  chosen  as  follows:  probability  of  detection  Pd  =0.8, 
probability  of  survival  Ps  =  1,  and  the  clutter  parameter  A c(z)  =  250^77  •  ^  ^as  n°ted  that  the 

assumed  observation  volume  is  uncertain  as  it  is  difficult  to  estimate  the  real  coverage  of  the  radar’s 
beam  suitable  for  filtering.  Qualitatively,  the  results  looks  a  good  fit  to  the  data,  but  for  the  filter  to  be 
confident  on  the  presence  of  the  object  it  takes  approximately  75  time  step  in  Figure  13a  and  150  time 
steps  in  Figure  13b.  In  between  steps  70  and  100  in  Figure  13a  the  filter  leads  the  estimation  to  the  wrong 
direction,  because  the  initial  estimate  obtained  after  a  relatively  short  target  present  around  step  70  was 
not  precise  enough  for  proper  propagation.  Visually  inspecting  the  data,  it  can  be  appreciated  that  the 
target  could  be  tracked  much  earlier.  For  the  second  data  set,  the  qualitative  results  are  noticeably  better 
with  the  filter  being  able  to  distinguish  strong  clutter  close  to  the  end  of  observation  (bottom  right  corner 
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(a)  Gaussian  mixture  (GM)-JoTT  filter  estimates  for 
the  Global  Precipitation  Measurement  (GPM)  satellite 
pass  (I). 


(b)  GM-JoTT  filter  estimates  for  the  GPM  satellite 
pass  (II). 


Figure  13:  Results  of  real-world  multiple-measurement  CAMRa  data  sets  using  the  standard  JoTT  filter. 


of  Fig.Figure  13b). 

Figure  14  shows  the  probability  of  detection  and  range  estimated  by  the  /3GM-JoTT  filter.  It  can  now 
be  appreciated  that  the  filter  does  manage  to  be  confident  to  have  an  estimate  much  earlier  compared 
to  the  standard  JoTT  filter  as  well  as  avoid  the  wrong  estimation  in  the  beginning  of  observations 
approximately  between  time  steps  70  and  100.  The  probability  of  detection  estimates  shown  in  the  upper 
plots  of  Figure  14  can  be  observed  to  correlate  with  the  presence  of  target-related  measurements  in  the 
range  plots.  It  can  also  be  noted  that  filter  continues  to  estimate  even  when  the  probability  of  detection 
Pd  drops  to  zero  and  no  target-generated  measurements  are  present. 

Finally,  the  results  using  the  r-JoTT  filter  are  shown  in  Figures  15a  and  15b.  To  evaluate  how  do  the 
probability  of  existence  q  and  the  probability  of  detection  are  related,  in  the  case  of  the  r-JoTT  filter,  the 
probability  of  survival  Ps  was  set  o  99.99%.  As  can  be  observed  from  the  figure,  the  estimated  number 
of  clutter  measurements  A©,  in  general,  is  below  zero  and  one.  The  explanation  for  this  is  the  result 
of  the  CFAR  processing.  CFAR  processor  produces  empty  detection  sets  for  the  observations  where  no 
appropriate  signal  strength  peaks  were  detected.  As  explained  in  [14],  the  majority  of  CFAR-processed 
radar  scans  are  empty,  thus  leading  the  A©  estimate  to  be  less  than  one. 

6.2  Vision-Based  Results 

6.2.1  Pre-processing  of  the  images 

An  important  pre-processing  step  to  perform  object  detection  and/or  tracking  is  the  reduction  of  noise  in 
the  images.  Bright  objects  such  as  stars  may  result  in  incorrect  detections  of  Near  Earth  Object  (NEO)s. 
As  explained  in  Section  5,  the  software  SExtractor  [3]  is  used  to  detect  stars,  but  also  returns  intermediate 
images  obtained  during  the  process  including  background  estimate  image,  background-subtracted  image, 
segmented  object  image  and  background-subtracted  image  with  objects  removed.  These  images  can  be 
seen  in  Figure  16,  the  negative  of  the  images  is  used  for  more  clarity. 

Adjusting  some  parameters  is  necessary  to  remove  most  of  the  stars  and  to  maintain  all  the  NEOs. 
DETECT_THRESH  is  an  important  parameter,  set  to  3.5.  It  indicates  the  value  at  which  the  pixels  are 
considered  objects.  In  this  case  the  threshold  is  3.5  times  the  standard  deviation  of  the  estimated 
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Figure  14:  Results  of  real-world  multiple-measurement  CAMRa  data  sets  using  the  /3GM-JoTT  filter. 


background.  DETECT  _M  IN  ARE  A  was  set  to  3  which  indicates  the  minimum  area  or  number  of  adjacent 
pixels  for  object  detection.  The  chosen  filter  was  a  3  x  3  top-hat  filter,  included  in  the  SExtractor  folder, 
tophat_2.5_3x3.conv.  The  filter  has  the  function  of  highlighting  interesting  pixels  and  reducing  the  noise 
level.  DEBLEND_NTHRESH  is  set  to  16  and  DEBLEND_MINCONT  to  0.00005.  Both  parameters  are  necessary 
to  determine  how  the  algorithm  splits  adjacent  pixels  belonging  to  different  objects,  recursively.  For 
background  noise  estimation  BACK_SIZE  is  set  to  24,  BACK_FILTERSIZE  to  5,  BACKPH0T0_THICK  to  24,  while 
other  parameters  are  set  to  their  default  values.  More  details  of  the  software  can  be  found  in  [3]  or  [9]. 


6.2.2  TLE  projection  over  the  images 

NEOs  are  seen  as  straight  lines  in  the  images.  The  length  of  the  line  depends  on  the  distance  from 
the  object  to  the  Earth  and  the  exposure  time  of  the  photograph.  The  TLEs  contain  information  of 
the  cataloged  debris  orbiting  the  Earth,  therefore  an  experiment  was  carried  out  to  project  the  objects 
presented  in  this  catalog  to  the  image  in  order  to  match  them  with  the  straight  lines  of  the  images.  The 
used  TLE  was  obtained  from  “http:/ /space-track. org”.  A  TLE  record  from  the  date  nearest  to  the  date 
of  the  observations  must  be  used  to  give  the  highest  accuracy  and  if  possible  from  the  same  day.  The 
PyEphem  library  [15]  from  Python  was  used  to  read  the  TLE,  and  predict  the  object  position  at  the 
same  time  the  photograph  was  taken.  PyEphem  internally  uses  the  SGP4  from  SGP  model,  and  provides 
accurate  astronomical  computations.  Given  a  date  and  location  on  the  Earth,  it  predicts  the  location  of 
satellites,  planets,  and  other  objects  in  different  coordinate  systems.  Highly  correlated  matches  between 
the  images  and  the  TLE  objects  are  shown  in  Figure  17.  Blue  stright  lines  correspond  to  the  projected 
debris  from  the  catalog  and  in  red  the  name  of  the  object  in  the  catalog  is  provided.  Note  that  the  blue 
lines  are  drawn  over  or  beside  the  debris  straight  line,  which  show  how  close  the  projected  and  the  real 
debris  are. 
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Figure  15:  Results  of  real-world  multiple- measurement  CAMRa  data  sets  using  the  r-JoTT  filter. 


7  Discussion  and  Future  Work 

During  this  period,  extensions  to  the  JoTT  filter  to  account  for  unknown  detection  statistics  were  devel¬ 
oped,  and  initial  work  on  detecting  orbiting  objects  from  visual  data  was  carried  out. 

Two  extensions  to  the  JoTT  filter  were  developed.  First,  to  account  for  an  unknown  and  possible  time 
varying  probability  of  detection,  the  /?GM-JoTT  was  developed.  Second,  to  further  account  for  unknown 
clutter  rate,  the  r-JoTT  filter  was  developed.  These  extensions  were  evaluated  using  simulated  data  and 
used  in  real  data  from  the  CAMRa  radar. 

Future  directions  related  to  the  JoTT  filter  are: 

•  Investigate  methods  to  account  for  unknown  birth  of  a  target. 

•  Implementation  and  evaluation  of  forward-backward  smoothing  [20]  and  the  extension  of  the  smooth¬ 
ing  algorithm  for  the  robust  versions  of  the  JoTT  filter. 

•  Comprehensive  approach  to  peak  (detection)  extraction  from  line  of  sight  measurements  of  CAMRa. 

•  Use  the  complete  radar  measurement,  i.e.  using  track  before  detect,  into  the  JoTT  filter  framework. 

In  order  to  advance  in  vision-based  RFS  tracking,  images  for  the  Falcon  project  are  waited.  In  the 
meanwhile  a  synthetic  database  will  be  build  based  on  the  OmniSSA  database  characteristics,  and  RFS 
algorithms  will  be  applied  on  them. 

7.1  Publications 

The  following  is  a  list  of  publication  related  to  the  Grant  No.  FA9550-15-1-0069. 
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(a)  Background  estimation  of  the  image 
(■ Ib(x,v )). 


(c)  Segmented  objects  from  the  image 
(■ Is(x,y )). 


(b)  Background-subtracted  image  y)). 


(d)  Background-subtracted  image  in  which 
stars  are  removed  (lR(x,y)). 


Figure  16:  Negative  of  images  obtained  during  the  process  of  object  detection  by  SExtractor. 
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(c)  Image  14,  Cosmos  2084.  (d)  Image  15,  Cosmos  2084. 


Figure  17:  Debris  match  between  the  TLE  catalog  and  the  image  straight  lines.  In  blue  the  projected 
debris  from  the  catalog,  in  red  the  name  of  the  object  in  the  catalog. 
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Acronyms 

/5GM  /3-Gaussian  Mixture. 


CA  Cell  Averaging. 

CAMRa  Chilbolton  Advanced  Meteorological  Radar. 
CFAR  Constant  False  Alarm  Rate. 

CPHD  Cardinalized  Probability  Hypothesis  Density. 

EKF  Extended  Kalman  Filter. 

FITS  Flexible  Image  Transport  System. 

GM  Gaussian  mixture. 

GPM  Global  Precipitation  Measurement. 

IOD  Initial  Orbit  Determination. 

JoTT  Joint  Target  Detection  and  Tracking. 

LEO  Low-Earth  Orbit. 

MAE  Mean  Absolute  Error. 

NEO  Near  Earth  Object. 

OmniSSA  Omnidirectional  Space  Situational  Awareness. 

PHD  Probability  Hypothesis  Density. 

r-JoTT  Robust  JoTT. 

RFS  Random  Finite  Set. 

RSOs  Resident  Space  Objects. 

SGP  Simplified  General  Pertrubations. 

SSN  Space  Surveillance  Network. 

TLE  Two-Line  Element. 

UKF  Unscented  Kalman  Filter. 

WCS  World  Coordinate  Systems. 
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Symbols 

x,  z  State  and  measurement  vector. 

F,  H  Transition  and  observation  matrices. 

Q,R  Transition  noise  and  measurement  noise  matrices, 
m,  P  Mean  and  covariance  of  a  Normal  distribution. 

X,  Z  State  and  measurement  set. 

Zi:k  Measurement  history  from  time  step  1  to  time  step  k. 

0  Clutter  RFS. 

q  Probability  of  existence  of  a  Bernoulli  RFS. 
s(x)  Spatial  distribution  of  a  Bernoulli  RFS. 

Pb,Ps,Pd  Birth,  survival  and  detection  probabilities. 

/b(x)  Spatial  distribution  of  the  born  target. 
fk\k-i{*k\xk-i)  Predictive  distribution  of  the  target  at  time  step  k  given  previous  state. 
fz( z|x)  Measurement  distribution  given  state. 
ft©( z)  Spatial  distribution  of  the  clutter. 

A©  Expected  number  of  clutter  measurements. 

G( A;  a,/3)  Gamma  distribution  with  parameters  a  and  {3. 

B(Pd]  a ,  b)  Beta  distribution  with  parameters  a  and  b. 

T(x)  Gamma  function. 

I(x,y)  Original  fits  image. 

I  block  Median  of  an  image  block. 

<j block  Standard  deviation  of  an  image  block. 

Ib(x,u)  Background  image. 

Background-subtracted  image. 

Is(x,y)  Segmented  image. 

th  Threshold  for  object  segmentation  on  Background-subtracted  image. 
Fc{x,y)  Gaussian  filter. 

Fth{x,y)  Top- hat  filter. 

Fb(x,y)  Block  filter. 

FMh{x,y)  Mexican-hat  filter. 
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